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ABSTRACT 


The thesis deals with the aeteeenen and tracking of dim 
point targets in infrared images. Research topics include 
image process modeling with adaptive two-dimensional Least 
Mean Square (LMS) and Recursive Least Squares (RLS) prediction 
Peeters. Target detection is performed by significance 
testing the prediction error residual. A pulse tracker is 
developed which may be adjusted to discriminate target 
dynamics. The methods are applicable to detection and 


tracking in other spectral bands. 
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I. INTRODUCTION 


A. PROBLEM DEFINITION 


In recent years considerable effort has been put into the 
improvement of infrared sensors. Increased sensor resolution 
was one of the more ardently pursued goals. The drive behind 
the desire for better resolution came from several areas. One 
of the driving forces was the military's need for infrared 
sensors that were capable of detecting heat signatured objects 
at long range. 

The specific concern in this thesis is the infrared 
detection of heat-signatured objects with sensors at altitudes 
of several hundred kilometers. Specifically, the thesis will 
be concerned with the detection and tracking of airborne 
targets such as medium range missiles and aircraft. The 
detection and tracking problem will be addressed from a 
generalized infrared sensor viewpoint. This is done with the 


desire to produce methods with the widest applicability. 


B. BACKGROUND 


The signature size and velocity of the targets have 
considerable bearing on the direction taken in solving this 
problem. Typical sensor geometric resolutions are on the 


order of .08 mrad. [Refs. 1 through 5] This produces a 


footprint on the earth's surface of approximately 48 m for a 
sensor at 600 km (1 Lrad=1 m at 1000 km). The type of targets 
considered in development of the problem have signatures 
smaller than the footprint of the sensor. This will be 
assumed to hold true even with the effects of point spreading 
and optical blur inherent in the sensors [Refs. 3 and 4]. 
Sensor sampling rates are typically on the order of 1 saneie 
per resolution cell, here .08 mrad [{Ref. 3]. The resultant 
intensity of a sample containing a target is due to the sum 
of target thermal intensity and the background intensity 
within the footprint [Ref. 2]. Thus the targets, with a 
modified intensity, in a discrete image are nominally about 
one pixel in dimension. 

The velocity of the targets is considerably larger than 
the drift rate due to sensor motion (e.g., two percent of the 
detector footprint) and jitter [Ref. 4]. If the image scan 
rate is typically 30 image frames per second, normal target 
movement between images will be less than two pixels per 
frame. A majority of the targets considered would actually 
move at less than one pixel per frame. 

Many of the infrared systems contain more than one focal 
plane assembly operating at different wavelengths. The 
resulting images can be dealt with as individual images each 
with its own statistical characterization. Ina single focal 


plane the scenes may be characterized by several 2-D random 


processes. The transition from one process to an adjacent 
process may be correlated or uncorrelated. The processes 
describe such items as terrain, and multiple layers of clouds 
and have a mean value dependent on the scene and the spectral 
band. (Ref. 4) Additionally the effects due to atmospheric 
Jitter, background scintillation, and electronic noise may be 
modeled by a combination of white and correlated eames ian 


newcse. 


C. CONCEPT AND CONTENT OF THE THESIS 

The sieesaculne sections have provided the background to 
develop a conceptual solution for the problem. It can be seen 
that ultimately the detection algorithm will need to detect 
pixel-size targets with intensities that may differ only 
slightly from the surroundings. Titec NG amr he. 71a ty, 
region of operation is seen to involve target motion of one 
to two pixels between successive image frames. 

In this thesis the solutions will be developed for images 
in a single focal plane. Their applicability will be general 
in nature and do not depend on the frequency band or the 
detailed sensor characteristics. Some ideas for merging the 
results from applying the algorithms in several focal planes 
simultaneously will be addressed in the final chapter. 

The next four chapters deal with the detection of the 


target and Chapter VI deals with the tracking. Initially a 


method will be proposed for modelling the background processes 
present in the image (terrain, clouds, etc.). In subsequent 
developments the background processes will be assumed to be 
homogeneous. Surveys of the literature show that frequently 
90% of the image may be characterized by a single random 
process. The modelling process will be performed with fixed 
length spatially adaptive two-dimensional prediction paves 
Adaptive filters used in the modelling process are the topics 
of Chapters III and IV. In these chapters the two-dimensional 
Least Mean Square (LMS) and Recursive Least Squares (RLS) 
filters are derived and tested. 

The detection process is the subject of Chapter V. Linear 
prediction produces an error residual process which is stored. 
This error process contains the unpredictable information in 
the image. In Chapter V a method based on statistical 
Significance testing will be used to separate the error 
residual into levels of predictability. The error residual 
present at the target location will be shown to be anomalous 
by comparison. This results in highlighting the targets after 
thresholding to achieve a low constant false alarm rate. 

In Chapter VI a point target tracker is developed and 
tested. The tracker is the based on a maximum a posteriori 
(MAP) estimate of the target position. In the overall system 
the detection algorithm passes its results to an autonomous 


tracking algorithm. Consideration of the logical interface 


between the detection algorithm and the tracking algorithm are 
addressed in Chapter VI. The precept throughout is the 
development of a set of algorithms that is reasonably generic. 
Algorithms developed extend readily to images other than 


infrared. 


II. MODEL 


A. GENERAL DESCRIPTION 


In order to pursue the objectives presented in the 
previous chapter it is necessary to establish a framework from 
which to work. Throughout this thesis the image to be 
processed is taken to be comprised of three components. The 
three components maybe viewed as a background generation 
process, a corrupting noise process, and a target generation 
process. Since target detection is the ultimate objective, 
it will be necessary to separate the three components. 

Background generation processes have been studied by many 
researchers; the result is a wide variety of models. If the 
image can be divided into statistically similar regions, then 
each region may be modeled on an individual basis. It will 
further be assumed that although the overall image may be 
nonstationary, the statistics in the separate regions are 
individually stationary. The region within the image is taken 
to be the result of a linear shift invariant operation on a 


wide sense stationary Gaussian distributed white noise source 
2 


of variance o-. [Ref. 6] A block diagram of the 2-D linear 


operation is shown in Figure 2-1. 


wom xtnmd 





Figure 2-1. Block Diagram 


The particular linear operation to be used is a_ two 


dimensional autoregression operation which may be written as: 


ime Smee, kK) x(n-!,m-kK) + w(n,m) 
& k 
Ce) ens 


(25) 
where $8 is some chosen region of support for the filter 
@eerficients. This model has been termed the White Noise 
Driven Representation (WNDR) and the transfer function may be 
viewed as an all-pole IIR filter with a constant numerator 
driven by a white noise source. The filter is said to be 
"Ccausual" (in the 2-D sense) if the computation of x(n,m) does 
not require knowledge of "future" values of the process. 
Otherwise it is said to be noncausal. The definition of 
causality is tied to the scan pattern used in processing the 
image. Figure 2-2 depicts the issue of causality for a row 


scan direction. 





Figure 2-2. Casuality Depiction 





If it is possible to determine the coefficients which describe 
the regional process, it then becomes possible to inverse 
filter the region. This will ideally leave only the driving 
noise process. 

Noise present in the image is assumed to be additive 
although not necessarily gaussian or white. If the noise is 
colored, the same procedure used to identify and eliminate the 
regional process would be utilized on the noise. The result, 
as before, is a white noise process. 

The most difficult of the three components of the image 
to characterize is the target process. If a specific model 
for the target is specified then the possibility exists that 


the overall system would be too restrictive. ine orders ro 


avoid this pitfall it is proposed that the only identifying 
features of the target are its pixel/sub-pixel size, 


opaqueness, and non-zero mean. 


B. PREDICTION ERROR FILTER 


It has been shown [Refs. 6 and 7], if the regional process 
is stationary and under specific conditions of process 
support, that an optimal inverse filter solution, in the MMSE 
sense, for the original process may be found. The drawbacks 
to the proposed solution are that the filter support with few 
exceptions will be of infinite extent. Further a priori 
knowledge of the 2-D correlation function or power spectrum 
and possibly the use of spectral factorization techniques are 
Beauiread to find the filter. Physical realization of the 
filter would, in any event, dictate truncation of the 
resultant filter support region prior to implementation. Once 
the filter support region has been truncated, it no longer is 
capable of matching the original regional autocorrelation 
function and may not be minimum phase and therefore not 
stable. [Refs. 8 and 9] 

The above considerations prompted the use of reduced order 
models for the regional processes and adaptive techniques for 
fitting the reduced order model. The support of the filter 
coefficients of the reduced order model is causal and has 


finite support as shown in Figure 2-3. 





Support 


The resultant reduced order model is 


x(n,m) = = ££ h(2,k) x(n-2,m-k) + v(n,m) 
Qk 
(2,k) € @ 


(2-2) 
where v(m,n) is a noise source, not necessarily white, that 
would be required to exactly reproduce the original regional 
process. 

The driving noise process may also be viewed as the error 
residual in the process produced by the inverse of the reduced 


order model operating on the original process as shown below. 


sae) 


e(n,m) = x(n,m) - = 5 h(2,k) x(n-2,m-k) 
Oak 
(laa oO 
(2-3) 

This representation is the prediction error filter (PEF) form 
(Ref. 9]. In practice if the support of the filter 
coefficients in the previous model is causal, and sufficiently 
large, the process complex spectral density function cs 
analytic in the neighborhood of the unit bicircle, and a 
minimum error variance fitting technique is employed to fit 
the PEF to the process; then the spectral density function 
can be closely matched. In this situation, the noise term 
e(m,n) is very nearly white and gaussian if the original 
process 1S gaussian. 

In practice the image background can only be assumed 
homogeneous (Gr. Ose stationary) in a small region. 
Consequently the model will have to be fitted to each local 
region or adapted during the image processing. Figure 2-4 
depicts a section of an image used for this processing. The 
point x(n,m) is the point at which the model of the process 
as presented in (2-3) is to be fitted. The prediction window 
represents the nonzero support of the summation term in (2-3) 
which will be referred to as the predictor. The prediction 


error filter is defined with a finite rectangular 
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Figure 2-4. 





Reduced Order Model Fitting 


support of size PxQ and for convenience of notation the 
coefficients are assigned as follows 


SUC) tate (Cyr — acre) 
a(@,k) = -h(@,k) oO < @ < P-1, 
0 < ae 


Q-1 


(2-4) 


The points underlying the prediction error filter mask may be 
concatenated into a column vector 


1 


earl, 1) 
aie 1) 


eet 1, m0) 
x(n, m=2®) 
oct 1 ) 


x(n,m-Q+1) 
x( n= leam—O-+s) 


Sn- P+ 1 moO! ) 


(2-5) 
The point ordering establishes the past (i.e., causality) with 
respect to the point x(n,m). From this vector the correlation 


matrix can be defined as 


(2-6) 
Since the region 1s assumed stationary, the correlation matrix 
is block Toeplitz. The optimal prediction coefficients are 
determined by solving the Normal Equations for linear 


prediction as shown below: 


(2-7) 


where i is a unit vector that is equal to one in the first row 
2 


e 


and zero elsewhere. The term o- represents the minimum error 
variance, in the MMSE sense, of the white residual error of 


the PEF using the optimal coefficients. In this thesis the 


3 


2-D Least Mean Square (LMS) algorithm is used to obtain an 
approximate solution to (2-7) in a local area. 

Under the assumption of local stationarity it is also 
possible to produce a set of prediction error filter 
coefficients for the region directly from the data. The 
method employed to produce that solution is to construct a 
family of equations of the form of (2-3) which represent the 
realization of the prediction error filter at each point in 
the region. These equations are then solved for the filter 
coefficients subject to the constraint that the sum of their 
error variances is minimized. The quantity to be minimized 
1s: 

3(x,a) = E(Z = [x(n,m) - X(n,m)}*) 
nm 
(n,m) € R 
(2-8) 
where x(n,m) is the prediction estimate (2-9) for the point 
at position (n,m) in region R, which has been previously 
defined as a region in which the process may be considered 


stationary, given by 


x(n,m) = = = h(£,k) x(n-£,m-k) 


Qk 
(2,k) € @ 
(£,k) # (0,0) 
(2-9) 
This procedure, known as "least squares," has been employed 
in many applications. In particular Ref. 10 has specific 
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bearing on this thesis. The recursive form of the least 
squares (RLS) technique has been used in the thesis. 

The final result of either the LMS or the RLS method is 
a "whitening" of the background process. The details of the 
actual processing will be covered in the remainder of the 


thesis. 


C. SIGNIFICANCE TESTING 


Under hypothesis testing a threshold is established which 
delineates the boundary between two hypothesises Ho and Hy: 
This can be seen to be a two sided testing procedure. The 
threshold is established by choosing the accceptable 
probability of rejecting Hy when Hp is true, P(type I error) 
and minimizing the probability of choosing Hp when H, Tot Lue, 
P(type II error). This method presupposes known distributions 
for the two hypotheses. It has been stated previously that 
a robust target detection system should minimize the 
Specifications on the target process. Hypothesis testing 
requires too much to be Known about the target. 

Significance testing is a one sided testing process with 
only one hypothesis Ho. refered to as the "null hypothesis." 
The procedure is to completely specify the probability 
PiscmibucLon Of Hy, then compute a test statistic from the 
event to be tested against Ho- The test statistic is a 


function of the event and is based solely on the event. The 


A ae. 


decision to accept or reject Ho after the statistic is 
compared to Ho is subjective, in that the establishment of 
boundary between the two decisions is at the descretion of the 
designer. Using the test statistic and the hypothesized 
probability distribpucion of Hj a region of values for the test 
statistic is determined where Hg will be rejected. This 
determines the probability of a type I error, which is known 
as the significance level of the test. The significance 
level may be based on desires to minimize the false alarm rate 
(time between type I errors), maintain a constant false alarm 
rate, etc. [Ref. 6] 

The significance testing used for target detection in this 
thesis involves a hypothesized gaussian distributed white 
noise background (Ho) - This noise is ideally the result of 
decorrelating the background as previously mentioned. The 


specifics of the testing will be covered in Chapter V. 
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IIIT. LMS FOR TWO-DIMENSIONAL PREDICTION 


A. THEORY 

In the previous chapter the use of a PEF with a reduced 
order predictor model was proposed as a method to decorrelate 
the background process. This, however, was contingent on 
Smeesing the predictor coefficients such that the error 
Variance at the output of the PEF waS minimized, see Figure 


el 


ee =I) | =pcOt 
© ) © 

= Gee ie Ge) 7) aiaG0, 27 
e ° © 


Figure 3-1. Prediction Error Filter 
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This chapter extends the traditional LMS algorithm to the 
problem of two-dimensional PEF. In this form the coefficients 
of the predictor are adaptively changed to minimize the mean 
squared error. (Refs. 11 and 12] The criterion to be 
minimized is 

€ = Efe(n,m;i)*] 

(sa 
The error term e(n,m;i), for the predictor mask of size PxQ 
with first quadrant support shown in Figure 3-1, is defined 
as 

P=170+1 

@ (i, m3 i) — oe ieee) => h(2,k:i) x(n-£,m-k) 
10) tere) 
(£,k) # (0,0) 

(3-2) 
where the index i represents the it” update of the filter. If 
the error is squared the criterion becomes 

, P=1 50s 
€ = E{x(n,m)"} - 2E{x(n,m) & eh (2, kei) x (ie —k)) 
t=0 k=0 
P=15,0-1 -P-17¢-1 
+E{ = ee De = h(g,h7;i) h(2,k:i) x(n-g,m-h) x(n-£,m-k) } 
g=0 h=0 2=0 k=0 
(g,h) # (0,0) 
(£,K) # (0,0) 

(3-3) 

The signal x(n,m) is assumed to be stationary and ergodic. 


It will be shown later that the stationarity assumption can 


be relaxed and still produce adequate results. The result of 
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(3-3) 18 a quadratic error surface that has a unique minimun. 
Differentiating (3-3) with respect to the predictor 


coefficients yields 


P= aol 
d€ = V(1) = -2E{(x(n,m) - = & h(g,hi7i) x(n-g,m-h) }- 
dh(2,k;i) g=0 h=0 
x(n-2,m-k) ) 
os 2 <¢ P-1 
O<k § Q-1 
(£,k) # (0,0) 


(3-4) 
This is the gradient of the error surface at update i mapped 
into the predictor coefficient space. The term x(n-£,m-k) 
represents the array of image points under the prediction 
mask; the gradient V(i) also is in array form. At the minimum 
point on the surface the gradient vanishes. The result of 
minimizing (3-4) and substituting into (3-3) produces the 
following expression for the minimum mean square error (MMSE): 

P=) sO=i 
é = E{x(n,m)*) - E(x(nm) EE h(£,k:i) x(n-2,m-k) } 


(3-5) 


B. THE STEEPEST DESCENT ALGORITHM (SDA) 
Central to the theme of the LMS algorithm is the method 


of steepest descent. This exploits the quadratic character 
of the error surface by iterative adjustment of the filter 


coefficients. The coefficients are adjusted by an amount 


= 


proportional and opposite in direction to the gradient at the 


current location on the error surface, that is 
h(@,k;itl) = h(2,k;i) + u(-V(i)) 
0o< 2 < P-l 
0O< k € Q-1 
(£,k) # (0,0) 
(3-6) 
Here, as in the one dimensional LMS [Ref. 11], the gradient 
is approximated by removing the expectation operator from (3- 
4) and recognizing the bracketed term as the prediction error, 


e(n,m;i), at update i. The form of the coefficient update (3- 


6) now becomes 


h(2,k;it1) = n(@,k;i) + 2u{e(n,m;i) x(n-2,m-k) } 
O< @ < P-l 
O< k € Q-1 
(£,k) # (0,0) 


(3-7) 
The term uw is a scalar known as the step size and controls the 


rate of adaption [Ref. 11]. 


C. ALGORITHM ANALYSIS 

In order to implement the LMS algorithm it necessary to 
use (3-2) and (3-7), and to choose an appropriate value for 
wu and an initial condition on the coefficient array. 

The LMS algorithm then attempts to track the error 
surface, under the assumption of stationarity, by adjusting 


the filter coefficients The resulee seaneeLronrsuarace seas 
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1s inherently changing even if the original image is 
stationary. 

In any analysis of IMS there is the problem of 
nonstationarity of the data over the entire image frame that 
needs to be considered. This problem has been dealt with from 
a variety of view points. In addition image processing has 
inherent problems caused by the overlapping of data blocks 
which was assumed. In the original derivation of the one 
dimensional LMS algorithm, independence of successive input 
data blocks is frequently assumed. This is not necessarily 
true for image data. This problem has been addressed in the 
literature under the assumption of specific dependencies. 
These range from independent input data to input data 
possessing considerable correlation. (Refs. 13-17] An 
expansion on the analysis 1s beyond the scope of this thesis. 
Experimental results indicate that Widrow's method of dealing 
with nonstationarity [{Ref. 11], as adapted by Bitmead [Ref. 
18] and Hsia [Ref. 19] and termed "The Normalized LMS (NLMS) ," 
has some general applicability in the processing of two- 
dimensional images. As originally presented the algorithm was 
based on the assumption of gaussian independent identically 
distributed (iid) data. In order to present an overview of 


this method, as adapted to two-dimensional processes, define 


2) 


<(n— em) 
x(n=2,m) 


® 


x (nN-P+1,m) 


x(n, m-O+1) 
x (n=1 ,m-Q+1) 


|x¢(n-P+1m-O+1) | 


(3-8) 
This vector consists of the inputs under the PxQ prediction 


mask at update i. Further define 


ill @) 
h(2,0) 


h(P-1,0) 


h(i) 


h(0,0-1) 
hii, Oa) 


h(P-1,Q-1) 


(3-9) 
as the prediction coefficients at update i concatenated into 
a column vector. Using the vector form of (3-7) it was 


proposed by Hadhound and Thomas [Ref. 12] that the actual 


coefficient adaption process in a image could be represented 


as 


Ze 


eee (ry ey (1) 

(om) 
The term y(i) is a column vector of random increments to the 
ideal time varying filter which would track the input data. 
Subtracting (3-10) from (3-7) in vector form one can define 

Gy) = h(i) 

(3-11) 
which represents the coefficient deviation from the ideal 
coefficients. The vector V(i) is called the coefficient error 
vector. The coefficient error update can be written as 
follows 

V(it1) = (I-2ux(i)x(i)')V(i) + 2ux(i)e°(i)-¥(4) 

(3-12) 
where e°(i) is the error resulting from prediction with the 
optimum time varying coefficients h°(i). The coefficient 
vector h°(i) represents the optimal coefficients at update i 
as determined by the solution of the normal equations. The 
ultimate desire is for the coefficient error term to approach 
zero aS 1 becomes large. If the input data is independent and 
stationary, this can be accomplished [Ref. 11] with a step 
size O<p<!/) , where Bene is the largest eigenvalue of R. 
However in the current form (3-12) it can be seen that this 
would not produce zero coefficient error as time progressed. 


The result is a contribution to the coefficient deviation due 


Zo 


to the approximation of the gradient (the first term of (3- 
12)) and a deviation due to the lag in tracking the ideal 
coefficient changes (the last terms in (3-12)) [Refs. 15 and 
16). 

In the "Normalized LMS (NLMS)" algorithm the step size is 
defined as 


aL 
‘ 4Tr(R) 





(355) 
where Tr(R) may be approximated by the (i) 2. The form 
given in (3-12) after substituting (3-13) has both heuristic 
and practical appeal. Both terms containing the input data 
under the mask have been normalized. This decreases the 
explicit dependence of the coefficient error on the change of 
data input power. The term a has been addressed in [Ref. 14] 
for the iid case and related to the ratio of the variance of 
the step change in (3-10) and the MMSE. This is of interest 
from a theoretical standpoint; however, in practice these 
factors are rarely available. In practice it has been 
observed that choosing a to satisfy O<a<2 for the case of 
independent data works well with the most rapid convergence 
of the MSE at a=1. In Ref. 14 an analysis was performed on 
correlated data for a two input predictor. The results 
extended to the "Normalized LMS" algorithm would indicate that 


the greater the correlation of the input data the smaller 
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should be the a@ term (typically a<l). After substitution of 
the step size (3-13) the coefficient update in vector form 


becomes 


a@e(n,m;1)x(1) 
e\(( ier JL))- 08)2\ (oe 2 er 
2Tr(R) 
(3-14) 
The coefficient error results in excessive MSE. A measure of 
this is the misadjustment ratio M defined as 
ae 
amin 
(3-15) 
where ¢€(©) is the MSE at steady state and € ie is the 
theoretical MMSE (Ref. 14]. 

Figure 3-2 shows a plot of the misadjustment ratio versus 
various values of the a term in the step size. The results 
presented are the average of simulations on synthetic images. 
The image process is modeled by a 2x2 quarter plane separable 
autoregressive filter. 


ene — 0.5 X(n—-1,m) — 0.4 X(mem-1) + 0.8°X(n=1,m-1) + w(n,m) 


(3-16) 
2 


Ww 


The results are for a range of o° from 0.6 to 4.0. It can be 
seen that this form of the step size produces a unique minimum 
M for a given process model. Considering (3-12) with Figure 


3-2, it can also be seen that aS a@ becomes small, the excess 
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MSE increases due to coefficient lag. Further,ge@as ia 


increases, so does the excess MSE due to gradient error. 


M (MISALIGNMENT) 
0.30 
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Figure 3-2. Misadjustment Versus a 


An example of the predictor coefficient deviation is shown 
in Figure 3-3. The image frame processed with the NLMS PEF 
was synthesized with (3-16) and o°=0.64. 

The initial coefficients were set by a least squares 
technique discussed in the implementation section. It can be 
seen for the processed Fon that the coefficients, while not 
the same as those for the original process, are stable and 
tracking. The row processed for Figure 3-3 is shown in Figure 


S=4. 
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Figure 3-3. Predictor Coefficient Plot 


PAY) 


0 


: 
: 


46 
ROW (PIXELS) 





Figure 3-4. Row of the Background 
Process 


The desired result is the decorrelation of the image 
process (1.e., a nearly white noise error residual); this can 
be seen in Figure 3-5. 

Figure 3-6 contains a plot of the normalized histogram of 


the error residual output of the NLMS Prediction Error Filter. 
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Figure 3-5. Error Residue Correlation 
Fonction 
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Figure 3-6. Error Residue Normalized 
Histogram Plot with Normal Density 
Function Superimposed 


29 


It can be seen that the error residual distribution has 
a gaussian appearance as expected with a gaussian driving 


process. 


D. LMS PEF IMPLEMENTATION 


One of the deficiencies in the LMS algorithm is the speed 
of convergence. The detection algorithm discussed ina later 
chapter depends on the ability of the PEF to decorrelate the 
background process. 

Another problem that is inherent in the processing of the 
image is that of directionality. In order to process the 
image a scan pattern and a method of handling boundary 
conditions must be chosen. Once these have been established, 
discontinuities and directional dependencies have been 
introduced into the process. Although directional dependence 
cannot be overcome, the effect of discontinuities and settling 
time can be decreased. 

The processing of images for this thesis was done using 
a PEF with a 3x3 rectangular support and a row raster scan 
pattern. Initial conditions were divided into two parts. The 
filter coefficients at the beginning of each frame were 
initialized with coefficients derived from a least Squares 
solution for a small area (9x9) pixels (see Figure 3-7) at the 
beginning of the frame. The initial coefficients could have 


been derived from other methods. However, the method used 
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Figure 3-7. Process Scan Diagram With 
Initialization Window 


should have a misadjustment less than or equal to that 
anticipated for the LMS algorithm. Initial coefficients for 
successive rows were steepest descent algorithm updates of the 
previous row initial coefficients. 

The remainder of the frame processing follows from the 
solution of (3-2) and the coefficient update computed using 
(3-7). The online computation of an estimate of the mis- 
adjustment ratio can be used as a measure of the algorithm 


efficiency. Anomalous jumps in the misadjustment can be used 
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to signal that the coefficients should be reinitialized using 


the least squares method. 


a2 


IV. RECURSIVE LEAST SQUARES FOR TWO-DIMENSIONAL 
PREDICTION 


A. THEORY 
ite wreast Squares @algerithm 1S a method used in the 
estimation of model parameters. (Refs. 6, 22, and 23] The 
objective is the estimation of the model parameters that will 
fit the model to the observed data. Mhe (Cri eer loner 
determining the goodness of fit is the minimum sum of the 
squared error between the model output and all observations. 
An extension of this method is the Recursive Least Squares 
(RLS) algorithm. This method is identical to the least 
squares algorithm but provides for updating the estimate of 
the model parameters as more data becomes available. im 
Chapter II the concept of fitting a reduced order model 
(prediction model) to an image was addressed. This model will 
be used in the RLS algorithm and has the form 
x(n,m) = > = hH(k,k) x(n-2,m-k) + e(n,m) 
Qk 
(2,k) € @ 
(£,k) # (0,0) 
(ered, 
The term x(n,m) is the observation and the terms x(n-2,m-k) 


are the past data. The terms h(£,k) are the unknown 


parameters that must be determined in order to fit the model 
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to the observation. The term e(n,m) is the error between the 
model and the observed image data at position (n,m). The 
least squares criterion is as follows 
idk Jyical 
J(h,xX;N,M) = & 5 e (h,X7nem) 
n=0 m=0 
(4-2) 
where J(h,x;N,M) is the cost function to be minimized by a 
appropriate choice of h(2,k;N,M). The error term in (4-2) is 
defined as 
e(h,xj;n,m) = x(n,m) - X(n,m) 
(4-3) 
The term X(n,m) is the prediction of the point x(n,m) defined 
as 
2a OS 1 
x(n,m) = 2 => h(’,k;N,M) x(n-£,m-k) 
£=0 k=0 
(£,k)#(0,0) 

Ga 
where the support of the predictor model is PxQ. The second 
indices (N,M) in the term h(£,k;N,M) are included to indicate 
that h(£,k;N,M) is held constant for a block of data 0<n<N- 
1 and O<m<M-1. Thus it is the sum of the errors with this 
choice for the coefficients that is minimized. 

The criterion as defined above in conjunction with the 


prediction model can be recast in the form of a PEF. Figure 


4-1 shows a 3x3 rectangular PEF as defined by (4-5). 
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Figure 4-1. Prediction Error Filter 





Peo = 1 
e(n,m) = x(n,m) - om => h(£,k;N,M) x(n-£,m-k) 
L=0 k=0 
(£,k)#(0,0) 
(4-5) 
The output of the PEF is the error term to be minimized as 
defined by the cost function. In order to simplify the 
notation the h and x variables have been dropped from the 
error term; however, the dependency is understood. As 
proposed in Chapter II, if the error term is minimized and the 


original process was a gaussian WNDR, then the error will be 


white and gaussian distributed. 
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At this point it is advantageous to define the previous 
expressions in vector form. Although the ultimate objective 
and results remain the same, directionality is inherent in the 
form to Gollow: The development to follow will be for a 
raster scan by row, but it may be generalized for any scan 
pattern for which recursive computability is assured. The 
Mapping function for a row scan (see Figure 4-1) is 

1 =n+ Nm 0O<n<N-1, O<n<N-1 
(4-6) 
where N represents the row length and M represents the column 
length of the data field defined on a rectangular lattice. 
By using this mapping function to replace the indices (n,m) 
and concatenating terms by rows into a column vector the 
following terms can be redefined. 


x (n-1,m) 
(On) 


x (n-P+1,m) 


x(n, m-O+F1) 
x (ns n-Ord) 


x (e+ moO) 


(A — ite) 
and 
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(4-8) 
ie wveCetor X(t), 1s thus comprised of the points under the 
prediction mask at update i=n+Nm. The vector h(j) is 
comprised of the prediction filter coefficients based on all 


data for O<i<j-1. It follows from this mapping that 


Xx(1) = x(n,m) and x(1) = X(n,m) 
(4-9) 
and 
&(i) = h'(j) x(4) O<isj-1 

(4-10) 

The cost function is then expressed as 
| eee, so Cl 
J(h,xij}) = 2d e(l)" = 24 lea ae hal esc) 
1=0 1=0 
(4-11) 


Where the term !“'*) is added to give the cost function a 
fading memory [Ref. 23] and thus a greater ability to adapt 


to nonstationarities in the image. Ina stationary image the 


ony 


exponential weighting factor would be set to unity. From this 
point on, the development of the RLS algorithm will parallel 
Haykin's one dimensional approach in Refs. 22 and 23. 

The cost function as presented is quadratic in the terms 
of h(j) and therefore possesses a unique minimum. To solve 
for the optimal coefficients h(j) in terms of the cost 


function, (4-11) is differentiated with respect to the filter 


coefficients. 
| J-1 j-(it1) Bees (ert |) 
d3(h,X7j) = -9 3} x(i) - =d h(j)x(1) 
dn (3) i=0 t= 


(4-12) 
Now as in Ref. 23 the following definitions will be made. 
The deterministic correlation matrix is defined as 


aa (1+ 1.) 


$(j) == x(i)x(i)! 
1=0 


(4-13) 
In the same vein the deterministic cross correlation vector 
(l1.e., the cross correlation of the point to be predicted with 
the data vector) is defined as 
ath gt digeal). 
8(j) = 22 Uy 
1=0 
(4-14) 
Using the last two definitions, setting the derivative of the 


cost function (4-12) equal to zero and rearranging yields 
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2@)7hya) 720) 
(4-15) 
This form has been termed in Refs. 22 and 23 as the 
deterministic normal equations. The solution of (4-15) 
produces the optimal coefficients h(j) in the least squares 


sense. 

h(j) = #1 (5) 803) 

(4-16) 

The objective now is to solve for h(j) recursively. This 

can be accomplished by forming the right hand side of (4-16) 

recursively. (Ref. 22] From (4-13) the correlation matrix 
can be expressed in the recursive form 

$(j) = d#(j-1) + x(5)x()! 
(4-17) 
Further, from (4-14), the cross correlation vector may be 


formed recursively as 
Cm eee (Iyer Cy PK yy 
(4-18) 
In forming the inverse of the correlation matrix the matrix 
inversion lemma [Ref. 22] is used. On application of the 


lemma, the inverse of the correlation matrix becomes 


x97! (5-1) x (5) x05) 1815-2) 
e105) = ate l-a) - 


1+ rxtxb (5) 64 (5-1)x(5) 


(4-19) 
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From (4-19) two terms are defined which will facilitate 
computations. These are 


P(j) = #1(35) 


(4-20) 
and 
x7) P(j-1) (3) 
kG) * | nn ae 
1+ A °xX (J) P(J-1)xX(9) 

(4=27) 
where K(j}) is called the "gain vector." The form of K(j) can 
be shown [Ref. 22] to reduce to 

ap) (a) ea) 

(4-22) 

It also follows from (4-19; 207-21 and 22) Venat 
P(j) = 7) pcj-ay - X 1 KS) x" (5) P(5-2) 
(4-23) 


The update of the filter coefficient vector h(j) can now be 


formed with the use of (4-16, 17, 21, and 22). 


Vg-1) x(9)] 


ND) a eal eee 
(4-24) 
Using this equation the prediction error is defined as 
e(j) = x(3) - bl (j-2)x(3) 
(4-25) 


This has been termed the a _ priori prediction error. Further 


define the a posteriori prediction error as 
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a(j) = x(3) - nl (3) x3) 
(4-26) 
It is the a prior prediction error that is used in the cost 
function in (4-2). Now using (4-25) in (4-24) the final form 
of the coefficient update equation is 
Reality + Ki) e(7) 
(4-27) 
This equation completes the necessary equations for the 
construction of the adaptively updated prediction error 
freer. 

A term that will be useful in the detection algorithm 
(covered in the next chapter) is the minimum value of the cost 
function. The cost function attains the minimum value if at 
each step in the recursion the filter coefficients satisfy the 
deterministic normal equations (4-15). In this case the cost 
function may be written as 

eye eet) 6 ( jac) 
(4-28) 
The weighted average of Juin(d) over the range of the 
iterations modified by the exponential weight factor provides 
an estimate of the error variance at step j. The effective 


memory length for averaging purposes is 


MEM = _A_ 
1-A 


(4-29) 
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B. ALGORITHM ANALYSIS 


The properties of the least squares method of parameter 
estimation have been addressed extensively in Ref. 23 and Ref. 
24. In Chapter II it was assumed that the background process 
could be modelled as a WNDR. It was further assumed that, if 
our reduced order model was sufficiently large and a minimum 
variance fitting technique was employed, the resultant error 
would be white. Under these assumptions the least squares 
estimate of h(2,k;M,N) has some important properties, namely 
that it is consistent and the best linear unbiased estimator. 
This has been shown for the one-dimensional least squares 
problem in [{Refs. 23, 24, and 25] and extends readily to the 
two-dimensional problem. Further it has been shown [Ref. 22 
and 23] for the one dimensional case that the mean square 
error will converge in approximately twice the filter length 
to the minimum mean square error. This was shown 
experimentally to be a good estimate for the two dimensional 
case as well. 

In addition, if the driving source for the model is 


gaussian, then the least squares estimate results in the 


maximum likelihood estimate. In order to demonstrate this 
some terms will need to be defined. Assume as before that 
the error residual is white and gaussian. The joint 


probability distribution for the error residual in (4-11), 


without the exponential weighting factor is 
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TO 2 
e 
(4-30) 
Inserting (4-25) for the term e(i) the likelihood function can 


be written as 


4-1 

Z (xi) - h'(3)x (i)? 
L Git 3) 0 ey = = 2 exp aa 450 
- (270.  i/ 20. 


(4-31) 
It can be seen that maximizing the likelihood of x(i), the 
observations, with respect to the parameter h to be estimated 
is the same as the least squares criteria. 

The algorithm has its faults as well as its good points. 
Without the exponential weighting factor the memory length is 
unlimited. This presents a problem in implementation on fixed 
word length devices. Normalization of the P matrix has been 
attempted to overcome this problem; however, this introduces 
a bias in the mean square error. The long memory also makes 


the algorithm less responsive to changes in the image as time 


progresses. This was the reason for the inclusion of the 
exponential weighting factor in the recursions. There is a 
drawback to the use of the weighting factor however. The 


smaller the factor (1.e., shorter the memory), the greater the 
error variance. Experimental results with synthetic images 


have shown this to be true. At a memory length of less than 
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approximately ten pixels the results were indistinguishable 


from those of LMS. 


C. RLS PEF IMPLEMENTATION 
The following is the complete RLS algorithm. Note that 


the initialization given here is the traditional one. The 
actual initialization techniques employed in this thesis are 
covered following the algorithm listing. The initial 
CONGLCIONS ON VP ye, sdamcsecae 


CO) oe 


1.0) 0 
Fo = 0 
(4-32) 
where 6 is a small number (approx. <<l). The exponential 


weight factor is set at the desired value in the range 0<i<1l. 
The update recursion for i>0 is then 

. 71 P(i-1)x(i) 
lS SS ee 
1 + lx! (i) P(i-2)x(4) 

(4-330) 
where K(i) is the gain vector. The inverse correlation matrix 
P is updated by 

P(i) = X!pc¢i-a) - X71 K(i)x' (i) P(i-2) 


(4-34) 
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The a_ priori prediction error is 


e(i) = x(a) - hl (i-2)x(i) 


(4-35) 
The filter coefficients are updated as follows 
iG) Ma sidiet) + K())e(5) 
(4-36) 
The posteriori prediction error update is 
a(i) = x(i) - h'(3)x(i) 
(4-37) 
mae final item to be computed is J,.. 
Seg jh ASeaeey (a= 1 ) + e(1)a(i) 
(4-38) 


It is from this last term that an estimate is formed of the 


error variance. This estimate is 
0° = 5,,,(i)/(n-8) 

(4-39) 
where n is the number of points processed and the n-8 term is 
to account for the loss of 8 degrees of freedom due to 
estimating the filter coefficients from the data set ([Ref. 
BAY 

The speed of convergence has been previously mentioned. 
While a convergence rate of twice the filter length may 
provide an acceptable rate of convergence for continuous data 


sets it presents a problem in image processing. The problem 


arises in three areas: 
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1) at the initial start up on each frame, 


2) at the beginning of each row or column in a row or 
column raster scan, and 


3) where a large change in background statistics occurs. 

In order to reduce the effects of the problem areas the 
least square algorithm was implemented in the form of a 
secondary processor. Logic flags were established in the RLS 
processing to signal the occurrence of one of the problems. 
Once one of the flags is raised the RLS algorithm is re- 
initialized with the least squares algorithn. The least 
Squares algorithm is applied to a 9x9 rectangular block of 
data centered on the flagged point in Figure 4-2. The 
algorithm produces an estimate of the filter coefficients, the 


P matrix, and J,.. which are passed to the main RLS processing 


n 
routine. 

It can be seen that this procedure can introduce a 
considerable burden in the processing. A suitable alternative 
is to eliminate the row initialization flag while still 
retaining the frame initialization flag. The initial estimate 
of P, h(j), and J|,, for each row or column could be obtained 
from the initial value of the previous row or column. The 


update is repeated for successive rows or columns in the same 


manner. 
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Figure 4-2. RLS PEF Initialization | 


An example of the experimental results follows using the 
two-dimensional RLS algorithm. The image processed was 
synthetically produced using a 2x2 autoregression mask. The 
actual process model was separable and defined as 
een = O.ox(N-t>m) — O0.4x(n,m-1l) + OVSx(n-1,m-1) F vi(n,m) 

(4=20) 
with the variance of the source aaa and output FoaA 
The overall mean of the image was removed prior to processing. 


The correlation function of the error residual resulting from 


processing the entire image is shown in Figure 4-3. 
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Figure 4-3. Error Residual Correlation 
Function 


Figure 4-4 is a histogram plot of the error residual with 
a gaussian Spoteclnoe superimposed for comparison. 

Figure 4-5 is a plot of the filter coefficients resulting 
from a RLS row raster scan of the same row as in Figure 4-3. 

This concludes the development and implementation of the 
two dimensional RLS algorithm. The RLS algorithm was found 
to be superior to the LMS algorithm in all areas except in the 
area of computations required. The performance will prove to 


be an overriding consideration in the target detection 


problem. 
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NORMAL DENSITY FUNCTION, N=3410 
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Figure 4-4. Histogram of Error Residue 
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Figure 4-5. Predictor Coefficient Plot 
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V. TARGET DETECTION 


A. GENERAL 

The previous three chapters have dealt with the theory and 
algorithms for modelling the background process. It was shown 
that the result of inverse filtering with an adaptive PEF is 
a nearly white error residual process. It is this error 
residual that is of particular interest in the detection 
algorithm. The residual represents the unpredictable portion 
of the image process. The question that arises is the effect 
of the resultant decorrelation on point targets present in the 
image. 

When a target is present in an image, it replaces the 
background pixel at that location. Therefore the target is 
opaque as far as the image process is concerned. Further the 
target, if present, is deterministic and its intensity remains 
constant over time. The probability density of the image 
versus the target may be viewed as in Figure 5-1. In this 
figure the intensity of the target is a fixed value T and its 
PDF is represented as an impulse; the mean of the PDF for the 
background is specified by m. It can be seen that only the 


Magnitude of the difference between the image mean and the 


Sil 


image => 


pixel magnitude 





Figure 5-1. Background PDF Versus 
Imbedded Target 


target intensity T is notable when compared to the background 
process. Therefore prior to processing with the adaptive PEF 
the background process mean is removed. This results in a new 
target T with an intensity of An. 

During the modelling phase the coefficients of the reduced 
order model are adapted to fit the image at each pixel 
location. The error in the fit is relatively small in the 
locations occupied by actual background pixels. This also 
holds true for any pixel value that closely resembles a 
background pixel; a target with a small Am is more difficult 
to distinguish from the background than one with a larger Am. 


The value of the error process at the target location is 


a2 


: : Pe a 
e (n,m) = T (n,m) - 2 > h(&,k)x(n-£,m-k) 
£=0 k=0 
(2,k)#(0,0) 
(5-1) 
The term e” (n,m) represents the error due to the presence of 
the target. This term may be rewritten as a combination of 
two error terms. The first term represents the error that 
would be result from predicting the background in the absence 
of a target. That is, 
P=iO= 1 
e(n,m) = x(n,m) - = & h(&,k)x(n-2,m-k) 
£=0 k=0 
(5-2) 
The other term is the difference between the target pixel and 
the missing background pixel. 
e7(n,m) = e (n,m) - e(n,m) = T (n,m) - x(n,m) 
(5-3) 
The resultant error process e"(n,m), separated in this manner, 
can be compared to the background error process. The process 
e(n,m), assuming complete decorrelation of the background, can 
be completely characterized by its variance. On the other 
hand the error process e7(n,m) possesses a mean which will be 
non-zero and is given by 
E{e7(n,m) }= Am = T" (n,m) 


(254) 
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The variance of error process e”" (n,m) is defined as 
Var{e (n,m)} = Ete” (n,m) °} - Ere” (n,m) }¢ 

(5-5) 
After using (5-1) through (5-4) in (5-5) and simplifying the 
variance of the target error process becomes! 

Var{e (n,m) } = o © (n,m) = o.*(n,m) 

(5-6) 
Intuitively this means that the ability to recognize the 
target is linked to the contrast between the background random 
process and the deterministic target. These two error 
processes present in a prediction error image may be viewed 
as in Figure 5-2. It can be seen that the magnitude of the 
mean difference between the target pixel and the anticipated 
background pixel determines the separation between the two 
density functions. It must be kept in mind that the density 
function is for the local area since the filter is spatially 
adaptive. Thus it is the difference between the local area 
mean and the target intensity that constitutes (Am) the 
difference between the two density functions. The factor is 
evident in an image with multiple background processes. The 
solution presented was for a homogenous (Single) background 


process. In an image with multiple background processes the 


lnote that as previously discussed x(n,m) and therefore 
e(n,m) has a mean of zero. 
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Figure 5-2. Background Error Process 
Versus Target Error Process 


removal of the image mean does not necessarily produce a local 
process PDF with a mean of zero. A separation between the 
local background process and the target process still results 
but it is now with reference to the local process mean. 

The form of the inverse filter as mentioned previously is 
a prediction error filter (PEF). This filter has a finite 
impulse response. It follows that the effect of a target on 
the error process is not limited to the pixel where the target 
is located. The PEF used in this thesis is rectangular with 
a finite impulse response of 3x3 pixels in dimension. If the 
background process is stationary within this region, the error 
process can be defined as follows. The error at each location 


within the impulse response support region is 


one, 


zZ Pa 
e"(nti,mt+j}) = & > a(l,kinti,mt+}) x(nt+i-2,m+j3-k) 
f=0 k=0 
(£,k)#(1,9) i 
+ a(i,msnti,m+jpe (n,m 
(357) 
where a(f£,k;nti,m+j) are the filter coefficients 


a(2,k;n+i,mt+}) = 1 (£,k) = (0,0) 
-~h(2,k) O<2<2, OK<k<2 


(5-8) 
The additional indices (nti,m+j) are to indicate that these 
coefficients are spatially variant. Using (5-3), (5-7), and 
(5-8) and concatenating the respective components into column 
vectors, the error process in the support region may be 


written as 


* 


e =ert ey;(n,m)a 

(S52) 
Now assume that the filter is at steady state within the 
region. Since the background has mean zero, e has mean zero. 
Then the mean of the error process at each location within the 

region in vector form follows from (5-4) and (5-9) 

Eje'} = Ama 

(5-10) 
Since the error process is uncorrelated, its covariance matrix 
K. is diagonal and follows from (5-6) and (5-8). The diagonal 


elements, k(i,j), of the covariance matrix are defined as 
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A) a a 
ata") Oo. +90, 1=)#0 
- (5-11) 
The purpose for the preceding derivation of the error 
process in the impulse region is twofold. First, it provides 
some additional information which could be used for methods 
of detection more oriented to hypothesis testing or 
significance testing of multiple pixel samples. Second, it 
shows that some significant error residual will occur in the 
immediate area of the target in the impulse response support 
region. The error is the result of the adaptive process that 
attempts to minimize the error variance. It seems reasonable 
that the greatest power in the error residual in the impulse 
response support region will be concentrated in the pixels 
which are the nearest neighbors of the target. This concept 


can be used in grouping suspected targets (results of the 


Significance test) and for track file reduction. 


B. SIGNIFICANCE LEVEL 


The concept of significance testing was discussed in 
Chapter II. Recall that a significance test is based on a 
null hypothesis Hy whose distribution is known. It is against 
this hypothesis that all other events are tested. The null 
hypothesis Hy to be used in this section is that of the error 


residual process. If the background process is gaussian then 


a) 


the error process is also gaussian. In the preceding section 
the error process was discussed. The error process due to the 
background was defined in (5-2) and the assumption was made 
that the process e(n,m) was completely uncorrelated. The 
probability density function for this process can thus be 


written as 


(alii) ) 2 + etn.m—,, 
(2m) /eq_ (n,m) 20, oth nye m) 


(5-12) 
The term o,(n,m) is indexed to reflect its dependence on the 
local statistics of the image. This density function can be 
used to define the significance level. 

Since the prediction error residual represents’ the 
unpredictable portion of the background it contains all of the 
pertinent information, without redundancy, necessary for 
making decisions. For this reason significance testing will 
be performed on the prediction error residual. The ordered 
set of all such residuals will be referred to as the 
prediction error image. 

It is desirable to set the significance level in such a 
way that little or no adjustment is needed over the entire 
image. To start with, a decision is made to reject Ho (i.e., 
to consider that a target may be present) whenever the 
magnitude of the error falls within a critical region. The 


lower boundary of this region is defined by the maximum 
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acceptable two tailed probability of a type I error, P (reject 
Ho given Ho is true) (see Figure 5-3). The computed 
probability of a type I error, based on the null hypothesis, 
will be called a. Associated with this probability is a level 


B which defines the critical region (see Figure 5-3). 


ao 0 B 
TD o< teritical reoion? 








Figure cae Diagram of Significance 
Testing Terminology 


Using these last two ideas the significance level a may 


be written as 


fo @) 
Qa = aa onn|- i ecnanl=|ae (nm 
(27) Oo, (n,m) 2 Oo, (n,m) 
B 


(5-13) 
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and the corresponding critical region is defined by 
jJe| > B 
(5-14) 
It can be seen that (5-13) and (5-14) depend on the local 
statistics. If a fixed threshold B (see Figure 5-3) is 
chosen, the probability a would fluctuate with the changes in 
variance o%(n,m) . This is undesirable. The probability a is 
equal to total area under the PDF within the critical region. 
The goal is to find a threshold which will produce a constant 
a independent of the local statistics. 
A solution is to map the error distribution into a 


normalized distribution with the mapping 


5u= (7m) 


a. (n,m) 

(5-15) 
This mapping results in a unit variance normal distribution 
in the variable z and the critical region is mapped into the 


region 


Al ees where A = —#— 
og, (n,m) 


(5-16) 


The probability a is determined by 
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2 
Ve) | —= aa exp 12° Jaz 
(2m) i/? 2 


(5-17) 
A depiction of the mapped error residual PDF is shown in 


Figure 5-4. 


HY) &X feriticalt region) 





Figure 5-4. Thresholding in the 
Normalized PDF Space 


Now observe that a choice of a threshold value of A=constant 
produces a constant probability a. The test (5-16) can thus 


be expressed, using (5-15) as 


2 
e{n,m > Xe = constant 
Oo, (n,m) 


(5-18) 
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From this it can be seen that a constant alarm rate (a) can 
be achieved with a constant threshold. The implication of 
this is that the statistic on the left hand side of (5-18) is 
independent of the variations in the background or error 
residual image. Ideally then, a single threshold would be 
Suitable for the entire image. This constant threshold may 
be found experimentally by adjusting the threshold in the 
prediction error images to produce on the average a specified 
percentage of false alarms (pixels which exceed the 
threshold). 

The preceding discussion defined the statistic to be 
tested and showed that the threshold may be set to any 
constant which produces the desired constant false alarm rate. 
The terms e(n,m) are computed using either the LMS or RLS 
method. In addition the term o° (n,m) needs to be computed. 
In the RLS algorithm this estimate was formed by (4-39). 
During the LMS process the estimate of o° (n,m) is computed 
from the difference in the unbiased spatial averages of the 
input x(n,m)¢ and the prediction g(n,m)¢. To compute this an 
unbiased estimate of the correlation matrix (R,) was used in 
conjunction with a column vector of the PEF coefficients as 
discussed earlier. The following equation defines the 
estimate of the error variance and is computed at each pixel 


location. 
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(5-19) 
This completes the description of the significance test. 
A block diagram of the process through the testing stage is 


shown in Figure 5-5. 


threshold 





Figure 5-5. Diagram of Detection Process 


The testing can be carried out at the same time the image is 
being whitened. This reduces the requirements for additional 
storage. It was mentioned earlier that the significance 
testing should be followed by a method for managing the 
resultant target list. This will not be addressed in this 
thesis. However a possible tracking algorithm is covered in 


the next chapter. 


C. SIMULATION RESULTS 


The simulations conducted using the processing algorithms 


developed up to this point were on images with homogenous 
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backgrounds. The image mean was computed and removed prior 
to processing; however, the local mean was not removed. The 
background of the image was generated from one of several 
models of the form 
2 2 
x(n,m) = 2% > h(2,k)x(n-2,m-k) + w(n,m) 

£=0 k=0 

(£,k)#(0,0) 

(5-20) 
where h(2,k) represents a fixed set of filter coefficients and 
w(n,m) iS a gaussian white noise source with a variety of 
variances. After preparing the image background, targets were 
placed at random locations within the image. The significance 
testing results with the various models and driving source 
noise were consistent throughout. The dimensionality of the 
PEFs used were 3x3 pixels; this was more than adequate for the 
2x2 model. It should be realized that in other cases the size 
of the PEF required to whiten a specific image may be larger 
or smaller than 3x3. 


The specific results presented in this chapter are for the 


background model in (5-20) with the coefficients h(1,0)=0.5, 
Z 


W 


h(0,1)=0.8, h(1,1)=-0.4 and o- =2.0. The imbedded targets 
were set to specific dB levels relative to the local variance 
in the 9x9 pixel region in which they were to be imbedded. 


This has been designated TBR (target to background ratio). 


Once the TBR was set, the mean of the 9x9 region was added to 
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the target. The target then replaced the pixel at the center 
of the 9x9 region. The purpose of this was to test the 
capability of the whitening and detection algorithms on 
targets that differed by very small amounts from their 
neighbors. Some simulations were ere for images with 
imbedded targets corrupted by additive white noise. The level 
of noise 1s specified as the BNR (background to noise ratio). 
The averaged results for multiple simulations at each 
parameter setting are shown in the Appendix. 

Plots of the results of processing identical images with 
the LMS algorithm and the RLS algorithm are shown in Figures 
5-6 and 5-7 respectively. These plots are for the conditions 
of no noise and a TBR of 5 @B with the threshold set to 11.4 
and ten targets imbedded in the image. Figures 5-8 and 5-9 
are the same images with imbedded targets of the same TBR as 
in Figures 5-6 and 5-7. However the images were corrupted 
with additive white noise. The BNR for both is 10 dB. These 
images were, respectively, processed with the LMS and RLS 
algorithms and the threshold was set to 11.4. 

Earlier in the chapter a model of the two error processes 
operating in a prediction error image was presented. It can 
be seen in Figure 5-2 that if the separation between the two 
processes is large, the threshold can be set to achieve a high 
probability of target detection while maintaining a _ low 


probability of false alarms. Refer to Tables 6 and 14 in the 
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Figure 5-6. Significance Test Results on 
an LMS Processed Image (No Noise) 
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Figure 5-7. Significance Test Results on 
an RLS Processed Image (No Noise) 
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Figure 5-8. Test Results 
Processed Image (BNR=10dB) 
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Figure 5-9. Test Results on ae RLS 
Processed Image (BNR=10dB) 


Appendix and recall the discussion (Chapter III) on the excess’ 
MSE resulting from use of the LMS algorithm. The RLS 
algorithm by its nature does not suffer from this type of 
excess MSE. Thus it would be expected that the whitened 
background process PDF would be broader for the LMS than for 
the RLS. This result can be seen in Table 6 and 14 by 


comparing the false alarm probabilities at which 100% 
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detection occurs. The false alarm probabilities essentially 
are a measure of the tail area in both tails of the background 
error process. It is seen from the tables that the broadened 
LMS error process adds an additional 9.25% to the tail area 
(false alarm probability). In another set of results, Tables 
3 and 11, the increased false alarm probability is seen to be 
approximately 6.5%. 

The effect of the increased background error variance can 
also be seen in Figures 5-6 and 5-7. In general from the 
tables it may be seen that at the same significance level the 
number of false alarms present remains approximately constant 
for both LMS and RLS. This result supports the proposed 
thresholding method and independence of the test statistic in 
(5-18). However the number of true targets found to be above 
the threshold is smaller for the LMS method than for the RLS 
method. 

The effect of the target error process in the impulse 
response support region can be seen in Figures 5-6 and 5-7. 
The scan direction for these figures was a row-wise raster 
pattern. The broadening of the error residual in the vicinity 
of the detected targets can be seen in the figures. It can 
be seen that the largest error residual occurs at the 
locations within one pixel of the true target location. 

It is possible to use the spread in the thresholds as 


an indicator of the spread in the target process. Bear in 
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mind that comparison between tables is subject to the effects 
Srecguaring the error in the™test statistic. Also™note in the 
tables the translation of the distribution of the thresholds 
is due to increases in the target magnitude. These points 
lend support to the assumed model. 

The preceding comparisons were the results of testing 
without added noise. It would be expected that adding noise 
would tend to increase the variance of the background error 
process while narrowing the target error process. The latter 
is due to the reduced ability of the adaptive filter to follow 
the changing background process. However in Tables 3 through 
5, for the LMS process with light noise (BNR=30 dB) the 
performance is actually better than with no noise. In 
adaptive gradient search algorithms, overlapping input data 
sets introduces correlation between the filter coefficient 
error vector and the data under the prediction mask. This 
creates an error bias in such algorithms as LMS [Ref. 13]. 
The addition of a small amount of white noise tends to reduce 
the correlation therefore decreasing MSE. This would decrease 
the width of the background error process and increase the 
variance of the target process. In Table 5 the beneficial 
effect is overcome by the increased variance of the target and 
background processes. Note that in Tables 11 through 13 for 
the RLS method, the background error process broadens as 


expected and as indicated by the increased false alarm 
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probability. In Table 12 some benefit is seen in light noise 
(due to target process shifting), but considerably less than 
occurred in the LMS case. A visual example of the effect of 
noise on the two methods is shown in Figures 5-8 and 5-9. 

From the plots and tables it can be seen that the RLS 
algorithm is superior to the LMS algorithm particularly in the 
area of dim targets and in the presence of noise. Of course 
the RLS algorithm is computationally more intensive than the 
LMS algorithm. The LMS computational advantage however may 
be offset by the additional computational burden of processing 
the increased number of false alarms. 

A point of particular interest is a region which may be 
called the "knee" of the significance testing process. This 
is the point at which the false alarm rate approximately 
doubles for a 10% increase in detection probability. This 
point appears in all of the tables and occurs at approximately 
80% detection probability for an LMS whitened image and 90% 
detection for RLS whitening. While these are not limits on 
the detection process they may be useful as figures of merit 


when comparing methods of whitening the image. 


D. CONCLUSIONS 
It has been proposed that the decorrelation of the 


background process results in a separation of the target 


process from the background process. The experimental data 


UZ 


supports this concept. Detection of targets at low false 
alarm rates requires an adaptive PEF that produces the 
smallest possible error variance. Thus it is not surprising 
in the comparison of the LMS and RLS methods that the RLS 
results are better. The false alarm rates using this 
detection method are relatively low for the detection rates 
achieved (see Appendix). This is due to the ability of the 
filter to adapt to local changes in the image. This in turn 
results in prediction coefficients which more closely model 
the local process and leads to a more complete decorrelation 
of the image by the PEF. 

The addition of white noise to images degrades the 
algorithm performance but not to the extent that might be 
anticipated. The noise broadens the background error process 
density function thus adding more false alarms for a given 
threshold setting. However the target error process variance 
is reduced and offsets some of the effect. Light noise 
(BNR=30 a@B) was found to increase the detection rate, 
particularly for LMS. Inevitably severe degradation occurs 
as the noise level increases above 10 GB. 

Significance testing showed some capability for detecting 
targets with TNR's near O GB. However the false alarm rate 
at that noise level may be higher than can be tolerated. The 
conclusion is that the methods used show considerable promise 


and capabilities in highlighting dim targets with TBR's of (1 


13 


to 2 dB) and in the presence of light to moderate noise 


levels. 
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VI. TARGET TRACKER 


A. GENERAL 

Throughout this thesis it has been assumed that the 
targets to be located are pixel or sub pixel in size. It has 
been shown that with the methods addressed in this thesis that 
target pixels can be identified. However this method also 
produces false alarms. The number of false alarms has been 
shown to be contingent on the significance level selected. 
If the significance level is set to achieve higher 
probabilities of target detection, the number of false alarms 
increase. Follow-on processing 1S required to reduce the 
number of false alarms and to track the targets spatially in 
successive images. 

The tracking method proposed in this chapter is based on 
the multi-dimensional pulse tracker proposed by Therrien in 
fete) 26). The tracker to be developed links the two 
objectives stated in the previous paragraph through the target 
motion. Three possibilities for the motion of the suspected 
target are no motion, uncorrelated, and correlated. Targets 
without motion will persist at the same location in the image 
through time. Uncorrelated motion is spatial displacement 


between images which is entirely random and unpredictable. 


HES) 


Correlated motion is random spatial displacement; however the 
target motion exhibits a trend in direction and ee Of 
movement. 

A tracker can be used as an observer, providing movement 
information that can be used to categorize the targets. This 
information and a_priori knowledge of the target dynamics can 
be used to reduce the target list. This is accomplished by 
comparing the observation to a specific profile of target 
dynamics. If the observation fits the profile, the target is 
retained in the target list; otherwise it is dropped. This 
thesis will not discuss the logic necessary to reduce the 
target list. The preceding was presented to indicate a 
possible means of using the tracker as a additional source of 
target information. Later in the chapter it will be shown 
that some target dynamic information can be incorporated into 
the tracker to help in tracking. 

Processing described in the previous chapters provides the 
location of possible targets and access to two types of image. 
These images are the actual observed image and the prediction 
error image. In these images the target has a deterministic 
shape which can be exploited in the construction of a tracker. 
It is assumed, as before, that the target in the observed 
image is approximately one pixel in size and has a non-zero 
magnitude. It was subsequently shown that the target, after 


whitening with the PEF, experiences point spreading due to the 
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filter's impulse response. An image containing a target with 
intensity profile s(n,m) can be modeled as 
r(n Mae) = s(nt+¢_,mt+¢ ;t) + w(n,m) 

(6-1) 
where r(n,m;t) is the observed intensity at a pixel (n,m) in 
an image at time t. The term s(nt+¢ ,m+¢ ;t) represents the 
target which should have occurred at (n,m) but Acer 
occurred at (n+¢ ,m+¢ ). The terms - and ¢, are the position 
error in the target location at time t. This can be viewed 
as a translation of the target from the expected location due 
to motion. The term w(n,m) is a additive noise term. The 
observed image in the absence of a target can be modeled as 

r(n,m;t) = w_ (n,m) + w(n,m) 

(6-2) 
where w_ (n,m) 1s colored noise (background or incompletely 
whitened background) and w(n,m) is additive white noise. The 
tracking problem is simplified if the term w_ (n,m) 1s near 
zero or if the target power is considerably larger than the 
colored noise component power. In either case the colored 
noise component can be ignored. For these reasons it was 
decided to perform the tracking in the prediction error image. 

During the background whitening stage the PEF coefficients 
are continually adjusted to produce the minimum error 


variance. This produces a near minimal phase filter and 
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results in most of the impulse response power occurring near 
the target location. With these considerations the target 
intensity s(n,m) in (6-1) was modeled as a spread gaussian 


pulse of the form 


C exp a —_ Oslin |<N=1 
2 


s(n,m) = O<|m|<M-1 





0 elsewhere 

(6-3) 
The term a controls the spread of the modeled signal. Fora 
given PEF impulse response support area, 0. is fixed. 
Experimentally the best value was found to correspond to be 
a radius around the target in which approximately 90% of the 
impulse response power occurs. For the adaptive 3x3 PEF this 
typically occurred at a radius of about one pixel. The scale 
factor C is used to match the signal magnitude and polarity. 
Polarity is positive if the target intensity is above the 
local image mean and negative if it is below. In order to 
match the model to the target signal, an estimate of the 
magnitude and polarity of the suspected targets are computed 
by the tracker when they are passed by the detection 


algorithm. 


¢although the true target intensity profile may not be 
Gaussian, this form is convenient. It will be seen later that 
this tracker is quite robust with respect to the detailed 
assumptions about target shape. 
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The detection algorithm identifies a possible target to 
the tracker. The tracker then extracts, from the prediction 
error image, a grouping of pixels centered on the target 
pixel. This will be referred to as the target tracking window 
and its center point will be called the centroid. A square 
tracking window is used with the method developed in this 
thesis. The window should be large enough to include most of 
the filter impulse response. Since two targets of equal 
magnitude falling completely in the tracker window would be 
treated as if there were only one, the window maximum 
dimension is based on the minimum distances between adjacent 
tracks. 

Once the target has been identified and the target 
location has been passed to the tracker, track processing 
becomes autonomous. The tracker will compute successive 
estimates of the target location in future images. With each 
estimate, the new tracking window centroid coordinates are 
computed and the tracking window is repositioned at the new 
centroid. This autonomous tracking method permits parallel 
processing for detection and tracking but may require a 
Separate prediction error image storage buffer for the 
tracker. The buffer is required if the processing speed of 
the whitening and detection processor is faster than that of 


the tracking algorithm. 


va, 


Additional logic would be needed to compare subsequent 
targets, identified by the detection processor, with existing 
tracks. Although not specifically addressed in this thesis 
some of the possible features of the logic would be: same 
target discrimination, split and merge of tracks, lost track 
identification, and so on. Ideally the complete tracking and 
logic package would produce a steady state number of false 
alarms tracks. This occurs when the number of new tracks 
started equals the number of old tracks terminated. A block 
diagram of the relation between the processing blocks in the 


detection and tracking system is shown in Figure (6-1). 





Figure 6-1. Detection and Tracking 
System Block Diagram 


The overall system, less connective logic, is a Pre-Whitening 
Unit (PWU), a Detection Processing Unit (DPU), a Tracking 
Window Control Unit (WCU), and the Tracking Unit (TU). The 


tracking window control unit and the tracking unit comprise 
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what has been referred to as the "tracker." The WCU controls 
the positioning of the tracking window based on the estimates 
in the target position provided by the tracking unit. The 
tracking window is then repositioned so that the centroid is 


at the estimated target location in successive images. 


B. TRACKER 

In the derivation of the tracker to follow the target 
pesition error (Co) in (6-1) is assumed to be a zero mean 
gaussian random process. Further it is assumed that the 
components of the error in the n and m direction are 
independent and that the noise in the models (6-1) and (6-2) 
is a zero-mean white gaussian random process. 

The origin of the coordinate system (n,m) is the centroid 
of the tracking window. The modeled target image (6-1) is 
located at an assumed position cn Se) which is to be 
estimated. A diagram of the tracking window and its 
PermenOoLogy 1S Shown in Figure 6-2. 

Initially an estimator for the position error of a single 
target pixel will be developed. This will then be extended 
to all the pixels in the tracking window. To begin with, some 
notation will be defined. Let the vector ¢ be defined as 

ae 


(6-4) 
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Figure 6-2. | Diagram of fTracking 
Terminology 


Further let z be the vector process of position errors in the 
target intensity at a pixel location in the tracking window, 


that is 


(6-5) 
Thus the vector z is comprised of the position error for a 
given pixel in the tracking window in each of k successive 


images. 
The vector r is Similarly defined as the sequence of image 


intensities at a given pixel location in the tracking window, 


that is 
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r(n,m;ti) 
(id; te ) 


| r(npmit*) 


(6-6) 
Finally the target vector s(¢) is a column vector containing 
the sequence of target intensity values present at the same 


location (n,m) in the tracking window, that is 


s(nt¢) mtg} tl) 
“ 2.4%) 
n m 


S(n+¢O -,mt+¢ 


k k -ak 
ig SES Bae ) 


(6-7) 
The index oe fou ei-l 2,35. k,1s che time index for the k 
image frames considered in estimating the position error. 
Each pixel intensity is a sample of the continuous Signal (6- 
3) within the tracking window. The observed target is assumed 
to be deterministic but subject to a random position error and 
additive white noise. The situation is depicted in Figure 
6= 2), 
The goal of the tracking algorithm is to predict the 


position error in the target in succeeding images. One method 
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of accomplishing this is to use MAP estimation. Given the 
statistical models for the position error process and assuming 
for the moment that the form of the target image process is 
known and given by (6-3) then it is possible to compute the 
a_ posteriori conditional density function p(z|r). From this 
density it is possible to determine the best estimate 2 in the 
MAP sense. It will be shown that future estimates of ¢ can 
be formed recursively from this estimate. [Ref. 26] Using 
Bayes rule the a posteriori probability density p(z|r) can be 


evaluated from 


(6-8) 
Thus using (6-1) the conditional probability density function 
Paes 2 aes 


p(x|z) = ——pfy—7pexpl- 2 (x-s(0)) RY (e-8 (0) 
(27) ane 2 


(6-9) 
where R. is the covariance matrix of the additive white noise 
process. The a priori probability density p(z) based on the 
assumptions stated earlier is 
+ 


a z'p oz) 
2 


Z) = 1 i 
p(z) (an) *72jR [eo 


(6-10) 
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where R_ is the covariance matrix of the position error in the 
n and m directions. Because of the assumed independence of 


¢ and ¢ this matrix can be partitioned as 
n m 


R" 


z 


0 


Ro =-— — 


R"™ 


| 
| 
a 
| 

(6-11) 
In order to find an MAP estimate we need to find the value of z 
that maximizes (6-8). Since the denominator is not a function 
of z it is sufficient to find z to maximize the numerator. 
If we substitute (6-9) and (6-10) into the numerator of (6- 
8), take the logarithm, and drop the constant term we are left 
with 


] 


TR 2) 


ecg) = $l (z-s()) RT (z-s(0)) + 2 


(6=12) 
To determine Zz) .. (the optimum estimate of z) the derivative 
of (6-12) is taken with respect to the position error vector z 


and the result is set equal to the zero vector. 
_9 m= as! (¢) =] — 
C = = R Eas =o z= 0 
oat Cfo) 78 , KS) , 2 


(6-13) 


Thus aa is defined implicitly by 


; 
5 ues ds =e 
a ie ie 0 Repe(Ssi(<) ) 


(6-14) 
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where the partial derivative term is a matrix whose components 
are [ds /oc. fle The signal present in any given image is 
subject to the position error in only that image. Therefore 


the partial derivative matrix has the form 


| 

Cn 2 | 

0 ds tt ) 0 ; 8) | 

cn | 

| 

eC) = | 
| 

0 as(c*) | 

| 


Cm ? 

0 ds(C_) 0 0 
oom 

0 : as(ck) 


(6=15) 
At this point it is beneficial to return to Equation (6- 
14). Now recall that the noise, w, was assumed to be white. 


Therefore the covariance matrix RL is a diagonal matrix with 


components equal to the noise power N\- Thus if we take 
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advantage of this and the partitioning of (6-15), the estimate 


for z can be written as 


2 (n,m) = nag ee (r-s(G)) 
: aoc |) Nie 
(6-16) 
and 
A (aya) = ny| aaa dA (r-s(£)) 
“L 90, | No 
(6-17) 


where Z (n,m) and Zz (n,m) are the position error estimates in 
the n and m directions, respectively, and are optimal in the 
MAP sense. 


Now define vectors a(n,m) and B(n,m) as 


a(n,m) = agi] 1 (r-s(£)) 
aC. No 
(6-18) 
and T 
B(n,m) = — 2  (r-s(£)) 
ac, | No 
(6-19) 


Then the two estimators can be written in the compact form 
z,(n,m) = R? a(n,m) 
(6=20) 


and 


(6-21) 


87 


It will be assumed for the moment that the covariance matrices 
R? and R® are known. Equations (6-20) and (6-21) then each 
represent a set of k scalar equations. The last equations of 


these sets have the forms 


p 
(2 (7p) ett 0 (np) 
: i=p-k+1 
(6-22) 
and 
¢ (n,m;p) = 5 RY (p-j) A(n, mi) 
i=p-k+1 
(6-23) 


th 


where a(n,m;1i) and B(n,m;i) are the 1°° components of a and £8 


defined as 


: 
a(n ine.) ae (r-s(C)) 
No oC 
(6-24) 
and 
B(n,mzi) = ca 6) (r-s(C)) 
No 
(6-25) 


These terms are the inner product of the gradient of the 
Signal multiplied by the difference of the observed and 
estimated images. These equations provide the means by which 
the target position estimate is updated. 

Consider the form of s(¢) Equation (6-7) in conjunction 


with the estimator Equations (6-22) and (6-23). It can be 
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seen that in order to compute the recursive estimates for the 
position error, k previous estimates are required. The upper 
limit in the summation represents the most recent image time 
1=p=0 representing the beginning of the track. It can be seen 
that the terms of (6-22) and (6-23) with indices i<p represent 
the track history prior to time i=p. These terms for a causal 
system are identically zero. At the beginning of the track 
the observed target location is passed to the tracker. Then 
the tracking window centroid is placed at the observed target 
location. By design the modeled target s(n+¢_,m+¢ ) is fixed 
at the center of the tracking window. Therefore the estimate 
of the position errors, ¢ and ¢ , at the beginning of the 
track is known and is incorporated in the centroid position. 
Thus all the essential terms are available recursively to 
compute the future estimates. 

The terms ¢ (n,m;p) and ¢ (n,m;p) pertain to a 
particular pixel within the tracking window. To compute the 
new position of the tracking window centroid, these terms are 
p+] 
n 


and APt! be the change in the 


m 


averaged. Specifically, let A 


centroid position. Then the change in the centroid position 


1s computed from 
p 
pt] _ 1 ee : 
A’ = 2 R_(p i). (2) 
1=p-k+1 
(6-26) 
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and 


apt 5 RT -i)B (i) 
m , (P avg 


1=p-k+1 
(6-27) 
where we have defined 
W-l1 W-l 
2 Z 
a > 2 a(n,m;1) 
04 1 = / 
Pet ae 
2 2 
(6-28) 
and 
W-1 W-l1 
2 2 
vyeh psi > 2 bin, tee 
a -— wate 
Pave! we n=-W+l m=-W+l 
2 2 
(6-29) 


The results of (6-26) and (6-27) represent the predicted error 
in the target location in the next image. The new centroid 


position is computed as 


centroid centroid n 
(6=3.05 
and 
pt+1 eee pt] 
Y entree Ine en BE A 
(G—= i 


The entire procedure is repeated in each succeeding image 
frame. In the absence of a priori knowledge of the future 


target position, the error in the new estimated target 
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location is assumed to be zero. It can be seen that the only 
information that needs to be retained to continue the 


recursion is the k past values of @ vg (t) and Cae 


C. TRACKER DYNAMICS 


The purpose of this section is to provide some insight 
into the dynamics of the tracker. Further a method of 
incorporating the target dynamics into the design of the 
tracker is shown. 

In the following analysis the motion of the target through 
the sequence of images will be assumed to be uncorrelated (see 
section A of this chapter). This type of motion is consistent 
with the problem of tracking small point targets with 
displacement between image frames of approximately one to two 
pixels. Since tracking is performed separately in the 
horizontal (n) and vertical (m) directions, we will only 
consider a single direction (horizontal) inthe analysis. All 
results shown also apply to tracking in the vertical 
direction. 

It follows from the assumption of uncorrelated target 
motion, that the error in the target position in the current 
image frame is independent of its position error in the 
previous frame. As a result, the position error correlation 


matrix R, is diagonal. The target position error estimates 


for the p+1°* image (6-26) and (6-27) then become 


ou 


7 
= 
L 
a) 


2 2 
Nee a » : ds(nt+o ,m+f_ )x"(x-s(C)) 
: NoWo  n=-W+1  m=-W+1 —o 
2 2 
(6-32) 
and 
W-1 W-1 
2 Pe 
 - oe z ds (n+ mt) E™(E-S(C)) 
NoW n=-W+1 m=-W+1 . 
2 2 


(6-33) 
where 0) and r” (0) are the values of the correlation 
function at lag zero. Since it was assumed that the position 
error process has zero mean, these terms are also the position 
error variances. Further note that the gradient is reduced 
to a scalar gradient of the modeled target image which is 
evaluated at each point in the tracking window. [In order to 
Simplify the notation in this section, the scalar gradient 
term in (6-32) will be written as V_ (n,m). 

With Equations (6-1) and (6-32) the mean of the position 


error estimate in the horizontal (n) direction is found to be 
n 
EA Pty = EO) Dae. V (n,m) E({r(n,m;p)-s(n,m) } 


Nowe n mm 
(n,m)ée{W,W} 


| 

RK 
n 3 
~_~ 
© 
—— 





ae? V_ (n,m) E{w(n,m;p) } 
Now’ n m 
(n,m) €{W,W} 


(6-34) 
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where w(n,m;p) is the additive zero mean white noise in the 


p'’ image frame. Further the variance of the position error 


estimate is 


n 2 
at (Cl eos = ag Ee V_ (n,m) (xr(n,m;p)-s(n,m;p) ) 
W n m 
. (n,m) €(W,W) 


ZS V (e,k) (r(2,k:p)-s(2,k:p))} 
‘i | 


0 
(£,k) €{W,W} 
: 2 
e(FAU) a,z 2 V¢nm? 
No Won m 
(n,m) €{W,W} 
: 2 
eyo PA) = (78) =e os 8 (n,m)? 


(n,m)e{W,W) 

(6-35) 

Note that if the tracker is to follow target position 

changes in the image, the variance of apt should be 

proportional to the variance of the change in the frame to 

frame position. Thus the tracker can be optimized to follow 

targets with a particular velocity by adjusting the term 
r’ (0). 

Figures 6-3 and 6-4 show plots of the target position 

error estimate versus the actual target position error. In 

Figure 6-3 the estimates are for target motion completely in 


the horizontal (n) direction. In Figure 6-4 the estimates are 
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Figure 6-3. Estimated Position Error 
Versus Actual Position Error; 


0.=2.0, L(G) a for .5 Pixel per 
Frame, EO) a leee for 1.0 Pixel per Frame 
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Figure 6-4. Estimated Position Error 
Versus Actual Position Error (Target 
Track 45. Degrees); o°=2.0, (OS o sexe 
.5 Pixel per Frame, Ea(O)=aoeeton ge. 0 
Pixel per Frame 


for target motion at 45 degrees away from the horizontal with 
a horizontal component of .5 or 1.0 pixels per image. In both 
plots it can be seen that the estimator response to the actual 
position error is divided into two regions. The first is a 
linear region centered at the origin and extending for about 
+1 pixel in true position error. The remainder is a region 
where the response is nonlinear. 

For these plots the term r’ (0) has been adjusted to place 


the known maximum displacement between frames in the linear 
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region of the estimator response. The region beyond this 
point is the nonlinear region. In this region the estimates 
of the target position error do not match the actual target 
position error. Thus, the result of setting the position 
error variance to match a specific target motion is that those 
with greater motion will tend to be rejected. 

It is apparent from the plots that only a small portion 
of the tracking window {-3,3} is contained in the linear 
region. This iS primarily a result of the gradient of the 
target model. Error contributions between the modeled target 
and the observed target (6-33) receive considerably less 
weight as the distance from the centroid increases. A 
possible means of extending the effective area of the 
estimator is through the use of a constant magnitude gradient. 


The form of the estimators for the "constant gradient" method 


is 
W-1 Wak 
atl = 4 Og -2sgn(n) x" (x-s(C)) 
NoW” n=-W+l m=-Wt+1 
2 2 
(6-36) 
and 
W- W-1 
iP 1 : : -2sgn(m)xr°(xr-s(¢)) 
NoW n=-W+i m=-Wtl 
2 2 
(6-37) 


where the gradient in (6-32) and (6-33) has been replaced by 
-2sgn(n) and -2sgn(m). Plots of the estimator response with 
the constant gradient are shown in Figures 6-5 and 6-6. The 
value of r. was decreased to produce estimators with the same 
power as those in the previous plots. The target motion in 
Figures 6-5 and 6-6 is the same as that in Figures 6-3 and 6- 


4 respectively. 
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Figure 6-5. Estimated Position .Error 
Versus Actual Position Error; o-=2.0, 
r"(0)=.1 for .5 Pixel per Frame, r'(0)=.2 
for 1.0 Pixel per Frame ; 
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Figure 6-6. Estimated Position Error 
Versus Actual Position Error (Target 
Track 45. Degrees); o°=2.0, fa) for 
.5 Pixel per Frame, r’=.2 for 1.0 Pixel 
per Frame 


It can be seen that using a constant gradient flattens the 
response of the tracker in the nonlinear region. The effect 
of maintaining the same power is to reduce in the magnitude 
of the estimate near the origin. However, this is offset by 
the increased magnitude of the estimates further from the 
origin. Thus, the effective response region has’ been 


increased virtually to the window limit (-3,3)}. 


98 


D. IMPLEMENTATION OF THE TRACKER 


After the images are processed by the target detection 
algorithm the resulting targets and prediction error image are 
passed to the tracker. The targets'-locations are placed in 
the tracker's "Track File" until retreived by the tracking 
window control unit. After retrieval of the target data the 
tracking window control unit accesses the prediction error 
image (see Figure 6-7). The data within the computed 7x7 
tracking window, centered on the designated target 


coordinates, is extracted. 


WINDOW 
CONTROL 
UNIT 


PREDICTION 
ERROR 
IMAGE 


TRACK FILE 


© 


Figure 6=7., 





Tracker and Trackfile 


At this point an estimate of the polarity and magnitude 


of the target is computed by 


a9 


3 3 
c=e(0,0)- 1 2 2 er 


46 n=-3 m=-3 
(n,m)#{ (0,0) (0,1) (1,0) } 

(6-38) 
where e(n,m) is the prediction error residue within the 
tracking window. The term e(0,0) is the target intensity, 
located at the tracking window centroid. LO MeONpuce =item. 
estimate the intensities of pixels considerably affected by 
the impulse response due to the target have been excluded. 
As mentioned previously these pixels are generally within one 
pixel of the target. 

After this computation, both the tracking window data 
array and the value of the scaling factor C are passed to the 
tracking unit. The tracking unit computes an estimate of the 
target position error using (6-32) and (6-33) or (6-36) and 
(6-37). The estimate is passed to the tracking window control 
unit which computes an updated tracking window centroid and 
repositions the tracking window. This new position and the 
estimated value of C are stored in the track file until a new 
prediction error image becomes available and the track file 
has been completely processed. 

A block diagram of the tracking unit is shown in Figure 
6-8. The summer computes the difference between the tracking 
Window array and the target model array and the output is 


multiplied by the reciprocal of the estimate for the noise 
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ARRAY MULTIPLIER 





Figure 6-8. Tracking Unit Diagram 


power. Two array multipliers are used to multiply this result 
by the target model gradient with respect to the n and m 
directions. For the general method represented by (6-26) and 
(6-27) the output of the array multiplier is averaged and the 
result is stored in a stack. The averaging is accomplished 
by the block designated as LPF. The oldest values of @ yg (t) 
and Boy, (4) are shifted out of the stack and discarded. The 
contents of the en and Csis buffers are vector multiplied by 
the k lags of the position error autocorrelation function. The 
result is the updated estimate of the position error. For 


uncorrelated target motion represented by (6-32) and (6-33) 


use of the stack is not required. 


uO 


Gradient and modeled target intensity values for locations 
within the tracking window are provided to both the summer and 
the array multiplier by the block labeled "Nonlinearity 8(¢)." 
This block is implemented as a look-up table. 

The system is initialized with the initial target 
coordinates, the previously estimated C term, and by setting 


@ vga? “Pee ee for all lags. 


FE. TEST RESULTS 


The following section is comprised of tracking runs at a 
variety of noise levels and for several different 
trajectories. All results were produced with the same 
uncorrelated trajectory model (i.e., only the zero-th lag of 
r, was not equal to zero). The results are presented in Table 
6-1 and in Figures 6-9 through 6-20. The table contains, for 
each figure, the signal to noise ratio, the estimated variance 
of the position error r (0), and the resulting sample error 
Variance (Var(n) and Var(m)) in the estimates. 

In the first trajectory (Figures 6-9 through 6-14) the 
target motion was .5 pixels per image frame in both the row 
and column direction. In the second trajectory (Figures 6-15 
through 6-20) the target was moving at approximately 1 pixel 
per frame. All figures in the simulations show motion 


through 20 image frames. This trajectory started at the lower 
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TABLE 6-1. TRACKING SIMULATION RESULTS 
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6-18(c) 


35 indicates that the pulse gradient method was used. 


4 indicates that the constant gradient method was used. 
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Figure 6-9. Pulse Gradient Model Track SNR=6.4dB 
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Figure 6-10. Pulse Gradient Model Track SNR=2.8dB 
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Figure 6-11. Pulse Gradient Track SNR=1.2dB 
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Figure 6-12. Constant Gradient Track SNR=6.4dB 
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Figure 6-13. Constant Gradient Track SNR=2.8dB 
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Figure 6-14. Constant Gradient Track SNR=1.2dB 
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Figure 6-15. Pulse Gradient Track SNR=6.8dB 
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Figure 6-16. Pulse Gradient Track SNR=2.8dB 
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Figure 6-17. Pulse Gradient Track SNR=1.2dB 
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Figure 6-18. Constant Gradient Track SNR=6.8dB 
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Figure 6-19. Constant Gradient Track SNR=2.8dB 
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Figure 6-20. Pulse Gradient Track SNR=1.2dB 
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left corner and moved to the upper right corner. In Figures 
6-9 through 6-14 a 180 degree direction reversal was made by 
the target, at the location (10,8), followed later by other 
abrupt changes in the trajectory. The figures show that the 
estimated trajectory follows the track well at all noise 
levels. The table further shows that the computed error 
variances, Var(n) and Var(m) are consistently less than one 
pixel. 

The estimates of the position error were computed uSing 
(6-3), (6-32), and (6-33) evaluated at the discrete points in 
the tracking window. This involved floating point operations 
and round off of the final position error estimate to the 
nearest integer value. As a result, quantization noise 
becomes a factor in the position error estimates. The 
quantization level separation in the case of an image is one 
pixel. Therefore assuming a uniform probability distribution 
for the floating point result, a quantization noise variance 
of .083 pixel¢ is expected. It can be seen for the higher 
Signal-to-noise level that the results are close to the 
quantization level. 

The noise rejection capability of the tracking algorithm 
using the pulse gradient form is seen to be high. In Figures 
6-9 through 6-11 and 6-15 through 6-17 the estimated tracks 
remain within one pixel of the true track. However, as might 


be expected, the constant gradient method is more susceptible 
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to noise and shows deviations from the track exceeding a pixel 
at even the higher signal-to-noise levels (see Figures 6-12 
through 6-14 and 6-18 through 6-20). However it may be 
required to use the constant gradient method if the target 
motion exceeds 2 pixels per image frame. When the target 
shift between frames exceeds the width of the window, the 
track will be lost. Simulations have shown that if only noise 
is contained in the target window the track will stagnate or 
Gere very little. The effective radius for the pulse 
gradient model was found to be approximately 1.50, where Oo. 
is the target model shape factor. This can be tied to the 
decreased magnitude of the gradient beyond this value. For 
the constant gradient method, the effective radius is 
essentially equal to the window dimension (W). However, 
because of this, the constant gradient method is more 
susceptible to clutter or targets crossing the boundaries of 


the tracking window. 


F. CONCLUSIONS 


This chapter was divided into several major areas. The 
first was a general description of the interface between the 
detection algorithm and the tracking algorithm. MThe other 
areas consisted of a development of the tracker equations, a 


dynamic analysis of the tracker, and simulation results. 


deals 


The results of these and additional simulations verified 
the tracker's expected performance. It was found through 
simulation that the gradient of the proposed target model 
limits the effective rate at which targets can move and still 
be tracked. This prompted consideration of the constant 
gradient model. In this model the modeled target intensity 
was left unchanged but the gradient was made a constant with 
appropriate sign. Use of the constant gradient was found to 
increase the track noise and decrease the resolution (i.e., 
minimum spacing between adjacent tracks) of the tracker. 

Either form of the tracker could be used as a target 
dynamics observer. An observer is a system which can produce 
estimates of parameters (e.g., velocity) that are not readily 
available. This observer information could be incorporated 
into the overall track file reduction process’ proposed 


earlier. This may be a profitable area for future research. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


A. GOALS 

The purpose of this thesis was to develop algorithms for 
the detection and tracking of dim point targets in infrared 
images. Images of this type include those produced by 
satellite-mounted infrared sensors. The targets are 
restricted to an approximate size of one image pixel and 
typically move at rates of approximately 1 to 2 pixels between 
successive image frames. A goal was also to develop these 
algorithms so they would be easy to implement and apply to 


other images in other spectral bands. 


B. SOLUTION 


The solution presented involves three steps. The first 
is to filter the image with an adaptive two-dimensional fixed 
size prediction error filter (PEF). This has the effect of 
whitening the image background and producing spikes at points, 
such as targets, different from the background. Both two- 
dimensional Least Mean Square (LMS) and 2-D Recursive Least 
Squares (RLS) filters were used. Both of these PEFs proved 
effective in whitening the background processes, however, the 


RLS method proved to be superior to the LMS method in the low 


1:9 


signal-to-noise ratio regime. This was shown to be important 
in the next step of processing. 

The second step involves significance testing the error 
residual produced by inverse filtering. This results in 
Suppression of the error residual produced by inverse 
filtering the background and highlighting the targets. The 
results presented showed that the LMS method could achieve a 
90% detection probability for 0 dB targets in the presence of 
light noise (30 dB SNR) with less than 9% false alarms. The 
RLS method produced comparable results for the same 
parameters. However, with moderate noise of (10 dB SNR), a 
2 dB target and a 90% detection probability, the LMS false 
alarm rate increased to 15% while that of the RLS increased 
only to 10%. The performance of the RLS detection algorithm 
was Significantly improved at target levels of approximately 
1 to 2 dB. In light noise this resulted in a 90% detection 
rate for an approximately 5% false alarm rate. Both detection 
methods were found to be easy to implement and relatively 
inexpensive. 

The final step taken in the problem solution was the 
development of a target tracker. This tracker receives 
initial target location information from the previous 
detection algorithm and tracks targets through successive 
frames of data. It was shown that a track can be maintained 


with less than one pixel error with signal-to-noise ratios 
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near OdB and for movement of 1 pixel per image frame. 
Additional simulations showed that track could also be 
maintained for 2 pixels per frame motion with approximately 
one pixel track error. The tracker was found to be robust 
with respect to the estimates of eee auieee parameters and 


to the general assumptions made in the development. 


C. RECOMMENDATIONS 

The previously mentioned processing involves images from 
only one focal plane. Often there are more than one focal 
plane, in different spectral bands, available in the infrared 
sensors. The significance of a target in one band will differ 
from that in the other bands. This will in turn result ina 
different PDF for the target error process. If the results 
of significance testing multiple images in different bands 
can be merged it may be possible to considerably reduce the 
false alarm rate. This is suggested as a possible area for 
further research. 

In order to implement the proposed method additional work 
needs to be done on the track file management logic. Some of 
the possible approaches and considerations were mentioned in 


Chapter VI Section A. 
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APPENDIX 


This appendix contains the results of the detection 
process using the LMS and RLS, algorithms. A detailed 
discussion is presented in Chapter V. 


TABLE 1. TEST RESULTS LMS PROCESSED DATA (NO NOISE) 
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pF is the ratio of the sum of non-target pixels above the 
threshold to the total number of pixels tested. 


Spd is the ratio of the sum of target pixels above the 
threshold to the total number of targets present. 
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TABLE 2. TEST RESULTS FOR LMS PROCESSED DATA 
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TABLE 3. TEST RESULTS LMS PROCESSED DATA 
(NO NOISE) 
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TABLE 4. TEST RESULTS LMS PROCESSED DATA 
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TABLE 5. TEST RESULTS LMS PROCESSED DATA 
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TABLE 6. TEST RESULTS LMS PROCESSED DATA 
(NO NOISE) 












r= == 


— | TBR (GB) way fe THRESHOLD 

re ee 
Pe | swe | ae | 
refs | oom | ao | a 
a 
en 
a 









127 


TABLE 7. TEST RESULTS LMS PROCESSED DATA 
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TABLE 8. TEST RESULTS LMS PROCESSED DATA 
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TABLE 9. TEST RESULTS RLS PROCESSED DATA (NO NOISE) 
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TABLE 10. TEST RESULTS RLS PROCESSED DATA 
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TABLE 11. TEST RESULTS RLS PROCESSED DATA 
(NO NOISE) 
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TABLE 12. TEST RESULTS RLS PROCESSED DATA 
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TABLE 13. TEST RESULTS RLS PROCESSED DATA 
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TABLE 14. TEST RESULTS RLS PROCESSED DATA (NO 
NOISE) 
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TABLE 15. TEST RESULTS RLS PROCESSED DATA 
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TABLE 16. TEST RESULTS RLS PROCESSED DATA 
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